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We propose a new approach towards determining the distribution of mechanical and acoustic wave 
energy in complex built-up structures. The technique interpolates between standard Statistical 
Energy Analysis (SEA) and full ray tracing containing both these methods as limiting case. By 
writing the flow of ray trajectories in terms of linear phase space operators, it is suggested here 
to reformulate ray-tracing algorithms in terms of boundary operators containing only short ray 
segments. SEA can now be identified as a low resolution ray tracing algorithm and typical SEA 
assumptions can be quantified in terms of the properties of the ray dynamics. The new technique 
presented here enhances the range of applicability of standard SEA considerably by systematically 
incorporating dynamical correlations wherever necessary. Some of the inefficiencies inherent in 
typical ray tracing methods can be avoided using only a limited amount of the geometrical ray 
information. The new dynamical theory - Dynamical Energy Analysis (DEA) - thus provides a 
universal approach towards determining wave energy distributions in complex structures. 



PACS numbers: 



I. INTRODUCTION 

Wave energy distributions in complex mechanical sys- 
tems can often be modelled well by using a thermody- 
namical approach. Lyon argued as early as 1967 23 that 
the flow of wave energy follows the gradient of the en- 
ergy density just like heat energy flows along the tem- 
perature gradient. To simplify the treatment, it is of- 
ten suggested to partition the full system into subsys- 
tems and to assume that each subsystem is internally 
in 'thermal' equilibrium. Interactions between directly 
coupled subsystems can then be described in terms of 
coupling constants determined by the properties of the 
wave dynamics at the boundaries of intersection alone. 
These ideas formed the basis of Statistical Energy Anal- 
ysis (SEA) which has since become an important tool in 
mechanical engineering and has been described in detail 
in text books such as by Lyon and DeJongSs, Keane and 
Priced and Craik^. 

A method similar in spirit but very different in its ap- 
plications is the so-called ray tracing technique. The wave 
intensity distribution at a specific point r is determined 
here by summing over contributions from ray paths start- 
ing at a a source and reaching the receiver point r and 
thus by the flow of ray trajectories. The method has 
found widespread applications in room acoustics^ and 
seismology^ as well as in determining radio wave field 
distributions in wireless communication 2 — and in com- 
puter imagining software 10 . A discussion of ray tracing 
algorithms used for analysing the energy distribution in 
vibrating plates can be found in 5 . 

Both methods - that is, SEA and ray tracing - are in 
fact complementary in many ways. Ray tracing can han- 
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die wave problems well, in which the effective number 
of reflections at walls or interfaces is relatively small. It 
gives estimates for the wave energy density with detailed 
spatial resolution and works for all types of geometries 
and interfaces. SEA can deal with complex structures 
carrying wave energy over many sub-elements including 
potentially a large number of reflection and scattering 
events albeit at the cost of reduced resolution. In addi- 
tion, the quality of SEA predictions may depend on how 
the subsystems are chosen as well as on the geometry of 
the subsystems itself, and a priory error bounds are often 
hard to obtain. 

Ray tracing and SEA have in common that they pre- 
dict mean values of the energy distribution and do not 
contain information about wave effects such as interfer- 
ence, diffraction or tunnelling giving rise to short scale 
modulations of the signal on the scale of a wave length. 
Both methods are thus expected to hold in the high fre- 
quency or small wave length limit where the small scale 
fluctuations in the wave solutions are often averaged out, 
for example, due to a finite resolution of the receiver. 

It will be shown here that SEA can be derived from 
a ray picture and is indeed a low resolution version of a 
ray tracing method. This makes it possible to introduce a 
new technique which interpolates between SEA and a full 
ray tracing analysis. The new method - Dynamical En- 
ergy Analysis (DEA) - keeps as much information about 
the underlying ray dynamics as necessary, benefiting at 
the same time from the simplicity of the SEA ansatz. 
DEA is thus an SEA type method in spirit but enhances 
the range of applicability of standard SEA considerably 
and makes it possible to give quantitative error bounds 
for an SEA treatment. 

The ideas as presented here have their origin in wave 
or quantum chaos theory in which short wave length 
approximations are combined with dynamical systems 
or chaos theory, see^ 2 - for an overview. Methods similar 



in spirit to the theory outlined in this paper have been 
discussed in the context of structural dynamics before. 
Heron 1 ^ modelled correlations between energy densities 
in subsystems which are not adjacent to each other in 
terms of direct and indirect contributions; the method 
does not take into account the actual ray dynamics 
and thus neglects long range dynamical correlations. 
Langley's 1 ^ Wave Intensity Analysis (WIA) treats the 
wave field within each sub-component as an (inhomo- 
geneous) superposition of plane waves thus introducing 
directionality which can propagate across coupling 
boundaries. The wave field is, however, assumed to 
be spatially homogeneous in each subsystem - an 
ad-hoc assumption which may often not be fulfilled. A 
ray-tracing treatment developed in a series of papers by 
Le Bot^ is probably closest to the approach presented 
here; by employing local power balance equations, a 
Green function for the mean energy flow is obtained 
and the full flow across subsystems is obtained via 
flux conditions. The approach differs in as far as we 
consider multi-reflection in terms of linear operators 
here directly and use a basis function representation of 
these operators leading to SEA-type of equations. 

The paper is structured as follows: in Sec. HH we will 
briefly review the ideas behind standard SEA. In Sec. 
Mil the ray tracing approximation will be derived start- 
ing from the Green function and using small wave length 
asymptotics. In Sec. HV] we will introduce the concept of 
phase space operators and their representation in terms 
of boundary basis functions. SEA emerges when restrict- 
ing the basis set to constant functions only. A specific 
example - coupled two plate system will be treated in 
Sec. El 

II. SEA REVISITED 

Starting point for an SEA treatment is the division 
of the whole system into subsystems; this will be done 
usually along natural boundaries, such as joints between 
plates or walls in a building. Energy is pumped into the 
system at localised or delocalised source points (such as 
the vibrations of a motor) and is distributed throughout 
the systems in terms of vibrational or acoustic energy 
in one form or another. The net power flow between 
subsystems is then given in the simple form 

Pi-j = WTHTIij , (1) 

V Hi rij J 

where Py is the power transmitted between subsystem 
j and i, uj is the (mean) frequency of the source, 
is the modal density of the (uncoupled) subsystem j, 
r/ij is a coupling constant and is the total vibra- 
tional energy in subsystem i. Allowing for a source term 
and dissipation and getting estimates for the coupling 
constant s 24 ' 27 ! 28 and the modal densities via Weyl's law, 
one obtains a linear systems of equations which is solved 
for the unknown energies Ej . SEA gives mean values for 
these energies in the same way as Weyl's law gives the 
mean density of eigenfrequencies. 



The validity of Eq. (jTJ) is based on a set of conditions. It 
is assumed that subsystems have no memory, that is, the 
coupling constants r/ij depend on the properties of sub- 
systems i and j alone. The eigenfunctions of the (uncou- 
pled) subsystems are furthermore expected to be locally 
described in terms of random Gaussian fields ("diffusive 
wave fields"). These key assumptions are expected to be 
valid only in the high frequency regime 2 ^ and for weakly 
coupled subsystems^ 7 .. 

By the nature of the technique, only relatively rough 
estimates for energy distributions can be obtained. Still, 
for high frequency noise sources, SEA or variants thereof 
are often the method of choice. 'Exact' solution tools 
such as finite element methods (FEM) become both too 
expensive computationally and unreliable, that is, small 
uncertainties in the systems may lead to very different 
outputs. One of the big challenges in mechanical engi- 
neering is the so-called mid-frequency problem - that is, 
handling the frequency range which is out of reach for 'ex- 
act' numerical methods but not yet in the high-frequency 
regime where SEA or ray methods are expected to work. 
SEA has been used as a starting point for penetrating 
the mid-frequency regime by employing hybrid-methods 
based on combining FEM and SEA treatments^ ' 19 ' 27 ' 30 . 

Connections between SEA and the properties of a ray 
dynamics associated with the wave equation has so far 
been made only indirectly. The statistical properties of 
wave systems with a chaotic classical ray dynamics have 
been shown to follow random matrix theory with wave 
functions behaving like random Gaussian waves^ 2 -. The 
basic SEA assumptions thus imply that the ray dynamics 
in each sub-element needs to be chaotic. This point of 
view has been stressed in the SEA literature by Weaver- 3 * 7 - 
and more recently in the context of determining the vari- 
ance of the wave output data in£ i 20 ' 21 . A detailed review 
discussing the connections between ray and wave chaos 
has been given by Tanner and S0ndergaard 3 ^. 



III. WAVE ENERGY DENSITY - FROM THE GREEN 
FUNCTION TO THE DIAGONAL APPROXIMATION 

A. The Green function 

We assume that the system as a whole is characterised 
by a linear wave operator H describing the overall wave 
dynamics, that is, the motion of all coupled subcompo- 
nents as well as damping and radiation. Different types 
of wave equations may be used in different parts of the 
system typically ranging from the Hclmholtz equation for 
thin membranes and acoustic radiation to the biharmonic 
equation for plate- like elements and to vector wave equa- 
tions describing in-plane modes in plates and bulk elas- 
ticity in isotropic or anisotropic media. We restrict the 
treatment here to stationary problems with continuous, 
monochromatic energy sources - generalising the results 
to the time domain with impulsive sources is straight for- 
ward. 

To simplify the notation, we will in the following as- 
sume that if is a scalar operator; treating bulk elastic- 
ity does not pose conceptual problems and follows for 



isotropic problems from Refi 3 - 3 - and for the anisotropic 
case, for example, from Ref. 34 . Note that both in SEA 
as in the new method - DEA- developed below, different 
wave modes such as pressure, shear or bending waves will 
be treated as different subsystems. 

The general problem of determining the response of a 
system to external forcing can then be reduced to solving 
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where the Green function G(r, ro, uj) represents the wave 
amplitude induced by a force Fo (of unit strength) acting 
continuously at a source point ro with driving frequency 
uj] p(r) denotes the material density. The wave energy 
density induced by the source is 
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(3) 



The bulk of the literature in acoustics and vibrational 
dynamics continues at that point by expanding the Green 
function in terms of eigenfunctions of either the full sys- 
tem or its sub-components. We propose to follow a dif- 
ferent route here by introducing a connection between 
the energy density and an underlying ray dynamics and 
expressing the Green function in terms of classical rays. 
A brief overview introducing the Eikonal approximation 
and the notation used for describing the ray dynamics is 
given in App. [A] 



B. Small wave length asymptotics of the Green function 



trajectory j at boundaries or material interfaces 2 ^ 16 ' 32 . 
Finally, A(g) contains geometric information and is of 
the form 



A (9) 
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where | • | = |det(-)| and the derivatives are taken in 
a local coordinate system r 1 - , r$ perpendicular to the 
trajectory at the initial and final point. The inverse 
of the Jacobian in Eq. ((6]) relates changes in the ini- 
tial momentum perpendicular to the trajectory, p$ , to 
changes in the final position r on the 'energy' mani- 
fold H = uj 2 = const. The phase index fij contains 
contributions from transmission/reflection coefficients at 
interfaces and from caustics, that is, singularities in the 
amplitude in Eq. ([6]). 

The representation ((4J has been considered in detail 
in quantum mechanics, see the books by Gutzwiller 1 , 
Stockmann 35 and Haake 13 . It is valid also for general 
wave equations in elasticity such as the biharmonic 3 
and the Navier-Cauchy equation 3 - 3 -; in the latter case, 
G becomes matrix valued. While the approximation 
is based on a small wave length expansion, the sum 
over trajectories in Eq. ((4]) often gives remarkably good 
results down to the mid and low frequency regime. Note 
that the summation in Eq. ((4j) is typically over infinitely 
many terms where the number of contributing rays 
increase (in general) exponentially with the length of 
the trajectories included. This gives rise to convergence 
issues, especially in the case of low or no damping, see 32 
and references therein. 



Using small wave length asymptotics, the Green func- 
tion G(r, ro,u>) can be written as sum over all classical 
rays going from ro to r for fixed H(r,p) — uj 2 , where H 
is the Hamilton function associated with the operator H, 
such as (|A2I) . and p is the momentum variable or wave 
number vector, see App. [AJ One obtains 1 - 

G{r,r ,uj)=C ^ a^jM-wt ^ 

j:r^r 

with prefactor 

r _ F Q n 1 

p(r )w(2 7 ri)( d + 1 )/ 2 ' 

where d is the space dimension. The action Sj(u>) is 
defined in Eq. (|A4[) . and is usually the dominant uj de- 
pendent term. The amplitudes Aj can be written in the 
form 3 *^ 

Aj = AfAfAf (5) 

containing contributions due to damping (d), conversion 
and transmission/reflection coefficients (c) and geomet- 
rical factors (g). The damping factor is typically of the 
form A)- = exp(— ctjLj) with aj, the damping rate and 
Lj, the geometric length of the trajectory. Furthermore, 

Ay corresponds to the product of reflection, transmis- 
sion or mode conversion amplitudes encountered by the 



The wave energy density, Eq. (|3|), can now be expressed 
as a double sum over classical trajectories, that is, 

e ro (r,uj)<x A 3 A r e l{s i- s j'- { ^-^' ) i'> (7) 

hi' -to— »r 

= p(r, r$,uj) + off-diagonal terms . 

The dominant contributions to the double sum arise from 
terms in which the phases cancel exactly; one thus splits 
the sum into a diagonal part 

p(r,r ,uj) = l^l 2 ( g ) 

containing only pairs with j = j' in Eq. (J7JI and an off- 
diagonal part containing the rest. The diagonal contri- 
bution gives a smooth background signal, which is here 
proportional to the energy density; the off-diagonal terms 
give rise to fluctuations on the scale of the wave length. 
The phases related to different trajectories are (largely) 
uncorrelated and the resulting net contributions to the 
off-diagonal part are in general small compared to the 
smooth part - especially when considering averaging over 
frequency intervals of a few wave numbers. (There are 
exceptions from this general rule; length correlations be- 
tween certain subsets of orbits can lead to important off- 
diagonal contributions. Coherent back-scattering or ac- 
tion correlations between periodic rays which have been 
identified to explain the universality of random matrix 
statistics are examples thereof; see 3 - 2 - for details). 



In what follows, we will focus on the diagonal part, 
that is, we will show that neglecting off-diagonal terms 
is equivalent to the standard ray tracing approximation. 
We will show furthermore that ray tracing can be written 
in terms of linear phase space operators and that SEA can 
be derived as an approximation of these operator. The 
connection between SEA and classical (thermodynami- 
cal) flow equations is thus put on sound foundations and 
the validity of the basic SEA assumptions as outlined in 
Sec. can be quantified. 



IV. PROPAGATION OF PHASE SPACE DENSITIES - 
FROM RAY TRACING TO SEA 

A. Phase space operators and probability densities 

We consider the situation of a source localised at a 
point ro emitting waves continuously at a fixed frequency 
lu. Standard ray tracing techniques estimate the wave 
energy at a receiver point r by determining the density 
of rays in r starting initially in ro (within the constraint 
H(ro,po) = lo 2 ) and reaching r after some unspecified 
time. This can be written in the form 

p(r,r ,Lj)=J drJdpJdX' (9) 

w(X',T)5(X-<p T (X'))p (X';u>) 

where X — (p, r) denotes a point ins phase space and the 
initial density 



Eq. (JHJ) can be simplified to 



p (X';uj) = S{r' - r )S(u 2 - H(X')), 



(10) 



is centred at the source point ro. Furthermore, X(r) = 
if T (X') is the phase space flow generated by equations of 
motion of the form (|A3j) with initial conditions X(Q) = 
X' and t is the time introduced in Eq. (|A3|) . It can be 
shown that Eq. ^ is equivalent to the diagonal approx- 
imation, Eq. ([8]), see App. [Bl 

The weight function w(X, r) contains damping and 
reflection/transmission coefficients and we assume here 
that w is multiplicative, that is, 

w(x, T 1 ) w (^ ri (X), t 2 ) = w{x, n + r 2 ), (n) 

which is fulfilled for (standard) absorption mechanism 
and reflection processes. Note, that the integral kernel 

C T (X,X r ) =w{X',t)8{X -ip T {X')) (12) 

is a linear operator - often called the Perron-Frobenius 
operator - which (after setting w = 1) may be interpreted 
as a propagator for the Liouville equation describing the 
time evolution of phase space densities^ 

p(X)) = {H(X),p(X)} 

(where {•, •} denotes the Poisson brackets) with solution 



p(r,r ,u) = 1 dr J dp w(p ,r ,r) 

6(r-rt(p>,r a ))5(Lu 2 -H(p',r )) 



(13) 



where <fil(X) — r(r) denotes the r-component of the 
flow vector. Eq. ([13")) is the starting point for a variety 
of ray tracing techniques popular in room acoustics^, 
seismology^ and optics such as illumination problems 
as well as for visualisation techniques in computer 
graphics^. 

While the basic equation (TT3"]) may seem 'obvious' from 
a ray geometrical point of view, we provide in Sec. IIII.BI 
and App. [B] a derivation from first principles starting 
from the wave equation. For references in a quantum 
context, se o 11 ' 29 . The connection between the ray tracing 
densities and the double sum over ray trajectories, Eq. 
([7]) , may form the basis for including " higher order" wave 
effects contained in the off-diagonal part. In what follows, 
we will stay within the diagonal approximations and use 
properties of the linearity of the phase space operator 
C to unveil the connection between ray tracing methods 
and SEA. 



B. Boundary maps and related operators 

We will for simplicity assume that the wave problem is 
confined to a finite domain with a well defined boundary; 
we may, for example, consider the vibrations of (cou- 
pled) plates of finite size or acoustics/elastic problems 
within bodies of finite volume. The long time limit of 
the dynamics is then best described in terms of bound- 
ary maps, that is, one records only successive reflections 
of a ray trajectory at the boundary. We introduce a co- 
ordinate system on the boundary, X s = (s,p s ), where s 
parameterises the boundary and p s denotes the momen- 
tum components tangential to the boundary at s; (X s is 
often referred to as Birkhoff coordinates). Phase space 
points X — (r,p) on the boundary are mapped onto X s 
by an invertible transformation B : X — > (X s ,ui) with 
H(X) = uj 2 . 

We now introduce two new operators: firstly, we define 
an operator Lb propagating a source distribution from 
the interior to the boundary, that is, 

C b {X s ,X')=w(X',t b ) cos0S(X s -B(cp TB (X'))) 

where X' is an arbitrary phase space point in the interior 
and tb {X') is the time it takes for a trajectory with initial 
condition X' to hit the boundary for the first time; the 
angle O(X') is taken between the normal to the boundary 
at the point s and the incoming ray velocity vector p, see 
Fig. Eh" Secondly, we introduce the boundary operator 

T(X S , A» = w(X' s )S(Xs - MK)), 

which is the Perron-Frobenius operator for the boundary 
map 



p(X,T) =C T [ Pa ] = J dX'5(X-^(X'))po(X'). ct> u (X' s )=B(<p TB (X')) with X' = B^iX'^uj). 



One can now write the stationary density in the 
interior, Eq. (fl3|) . in terms of the boundary opera- 
tors introduced above. Firstly the initial density fTO]) 
is mapped onto the boundary, that is, po(X s ,u>) = 
J dXCsiXs, X)po(X,lu). The stationary density on the 
boundary induced by the source po(X s ,u>) is then 

oo 

p{u) = T n {u) p (u) = (1 - r( W ))-^(w) . (14) 

n=l 

where T n contains trajectories undergoing n reflections 
at boundary. The resulting density distribution on the 
boundary, p(X s ,u>)), can now be mapped back into the 
interior using C^ 1 and one obtains the density (fT3|) after 
projecting down onto coordinate space, that is, 



p(r, r ,v) 



J dpdX s C B \X,X s )p{X s ,Lo) 



(15) 



The long term dynamics is thus contained in the op- 
erator (1 — T) _1 and standard properties of the Perron- 
Frobenius operators ensure that the sum over n in Eq. 
([T4| converges for non-vanishing dissipation. Note, that 
for w(X) = 1, T has a largest eigenvalue 1 and the ex- 
pression in Eq. (fT4]) is singular. That is, in the case of 
no losses due to absorption or radiation, a source con- 
tinuously emitting energy into the system will lead to a 
diverging energy density distribution in the large time 
limit. The eigenfunction of T (and C) corresponding to 
the eigenvalue 1 is the constant function; that is, in equi- 
librium the energy is equally distributed over the full 
phase spaced. 

To evaluate (1 — T) _1 it is convenient to express the 
operator T in a suitable set of basis functions defined on 
the boundary. Depending on the topology of the bound- 
ary, complete function sets such a Fourier basis for two 
dimensional plates or spherical harmonics for bodies in 
three dimensions may be chosen. Denoting the orthonor- 
mal basis {^o, ^2, • ■ ■}) we obtain 

T nm = J dX s dX' s ^* n (X s )T{X s ,X' s ;u;)^ m (X' s ) (16) 



dxi 



,(X' s ))w(X' s )^ m (X' s ). 



The treatment is reminiscent of the Fourier-mode approx- 
imation in the wave intensity analysis (WIA)i&; note, 
however, that the basis functions cover both momentum 
and position space and can thus resolve space inhomo- 
geneities unlike WIA. If the boundary map 4> U {X S ) is not 
known or hard to obtain, it is often convenient to write 
the operator in terms of trajectories with fixed start and 
end point s' and s; one obtains 



T„ 



dsds' 



= I dsds' 
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\ds/d P ' s 
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dsds' 



-K(X s )w(X' s )^ m (X' s ) (17) 
K(X s )w(X' s )^ m (X' s ) 



with X s = (s,p s (s,s')) and X' s — (s',p' s (s, s')) and S 
is the action introduced in Eq. (|A4[) . The representa- 
tion, Eq. ([17) , is advantageous for homogeneous prob- 
lems where the ray trajectory connecting the points s' 



and s is a straight line, see the examples discussed in 
Sec. El 



C. Subsystems 

In many applications, it is natural to split the full 
system into subsystems and to consider the dynamics 
within each subsystem separately. Coupling between 
sub-elements can then be treated as losses in one sub- 
system and source terms in the other. Typical subsys- 
tem boundaries are surfaces of reflection/transmission 
due to sudden changes in the material parameters or 
local boundary conditions due to for example bends in 
plates. Also, weakly connected sub-domains such as two 
regions connected through small openings may be consid- 
ered as separate subsystems. We denote the subsystems 
{Pi, . . . P/v} and describe the full dynamics in terms of 
the subsystem boundary operators T lJ ] flow from Pj to 
Pi is possible only if the two subsystems are connected 
and one obtains 



(18) 



where <p l J is the boundary map in subsystem j mapped 
onto the boundary of the adjacent subsystem i and X l s 
are the coordinates of subsystem i, see Fig. [T}d. (Note, 
that subsystems exchanging wave energy are necessar- 
ily connected through a common boundary here). The 
weight w lJ contains, among other factors, reflection and 
transmission coefficients characterising the coupling at 
the interface between Pj and Pj. 

A basis function representation of the full operator T 
as suggested in Eq. (|16p is now written in terms of sub- 
system boundary basis functions with 



dXtdXi ^ n (Xi) TV(Xt,Xi) «!';„;. Ym. (19) 



The equilibrium distribution on the boundaries of the 
subsystems is then obtained by solving the systems of 
equations (fT4)) 



(1 - T)p = po 



(20) 



Here, T is the full operator including all subsystems and 
the equation is solved for the unknown energy densities 



P={p\ 



where p l {n) denotes the (Fourier) coef- 



ficients of the density on the boundary of subsystem i. 
Equations similar to (f2"0"]) have been considered by Craik 5 
in the context of SEA. Note, that for a source localised in 
subsystem j, one obtains p l ^ only if Pi has a boundary 
in common with subsystem Pj. 



D. From ray tracing to SEA 

Up to now, the various representations given in Sec. 
IIVI are all equivalent and correspond to a description of 
the wave dynamics in terms of the ray tracing ansatz §§§ . 
Traditional ray tracing based on sampling ray solutions 



FIG. 1. Coordinates used for the boundary maps: a) in case of a single sub-system; b) at an intersection between two 
sub-systems 



over the available phase space is rather inefficient, how- 
ever. Convergence tends to be fairly slow, especially if ab- 
sorption is low and long paths including multi-reflections 
need to be taken into account. Finding all the possible 
rays which connect a fixed source and receiver point is a 
computationally expensive boundary value problem and 
typically only a small sample of all the trajectories cal- 
culated are actually needed in determining the local en- 
ergy density. In addition, the number of rays connecting 
source and receiver grows quickly (often exponentially) 
with the length of the ray trajectories setting fairly tight 
numerical bounds on the number of reflections one can 
take into account - a severe limitation in the low damping 
regime. 

These problems are common for ray summation 
methods^. They can be overcome by describing the 
dynamics in terms of boundary operators and boundary 
functions ^ n as outlined above. While the representa- 
tions are equivalent when employing the full set of basis 
functions (leading to infinite dimensional operators T), 
this is, of course, not the case for finite dimensional ap- 
proximations. When considering the solutions of Eqs. 
(fT4|) or (|20|) . one is in general interested in smooth ap- 
proximations of the energy density obtained from the 
classical flow. The resolution required is naturally lim- 
ited by the wave length of the underlying wave equation, 
but in many application a much coarser resolution will be 
sufficient. Convergence for obtaining such coarse grained 
energy density distributions is in general fast when in- 
creasing the dimension of the operators involved and of- 
ten only a very small number of basis functions (of the 
order < 10 per subsystem and momentum and position 
coordinate) are necessary. In addition, only short ray- 
segments are needed to evaluate operators of individual 
subsystems as multi-reflections are included explicitly in 
the sum (|T3|) . 

An SEA treatment emerges when approximating the 
individual operators T tJ in terms of the lowest order ba- 
sis function (or Fourier mode), that is, the constant func- 
tion \&q = (Ag)^ 1 / 2 with A B , the area of the bound- 



ary of Pj. The matrix T y is then one-dimensional and 
gives the mean transmission rate from subsystem Pj to 
Pi. It is thus equivalent to the coupling loss factor rjij 
used in standard SEA equations (JTJ) - The resulting full 
iV-dimensional T matrix (with N, the number of subsys- 
tems) yields a set of SEA equations using the relation 
PO)) (after mapping the boundary densities back into the 
interior with the help of local operators C B ). Note, that 
the terms Ei/rii ~ p 1 in Eq. (|TJ) are in fact mean energy 
densities as the mean density of eigenvalues is to leading 
order rii cx Ai with Ai the area/volume of subsystem Pj 
following Weyl's law^. 

The matrix T can in this approximation be interpreted 
as a transition matrix of an N dimensional Markov 
chain; SEA is thus in fact a Markov approximation of a 
deterministic dynamics. Similar approaches have been 
taken in dynamical systems theory over the last decades 
leading to a stochastic interpretation of chaotic dynami- 
cal systems in terms of a thermodynamical formalism 7 . 
A Markov or SEA approximation is thus justified if 
the ray dynamics within each subsystem is sufficiently 
chaotic such that a trajectory entering subsystem j 
'forgets' everything about its past history before exciting 
Pj again. In other words, correlations within the dy- 
namics must decay fast on the timescales of the staying 
time fj, that is, on the time scale it takes for a typical 
ray to leave Pj either by being transmitted to another 
subsystem P t or by being lost due to absorption. The 
dynamics must indeed equilibrate on the time scale fj. 
This condition will often be fulfilled if the subsystems' 
boundaries are sufficiently irregular, the subsystems 
are dynamically well separated and absorption and 
dissipation is small - conditions typically cited in an 
SEA context. In this case, SEA is an extremely efficient 
method compared to standard ray tracing techniques. 
However, for subsystems with regular features, such as 
rectangular cavities or corridor- like elements, incoming 
rays are directly channelled into outgoing rays thus 
violating the equilibration hypothesis and introducing 
memory effects. Likewise, strong damping may lead 



to significant decay of the signal before reaching the 
exit channel introducing geometric (system dependent) 
effects - that is, the distance between input and output 
channel becomes relevant. 

These features can all be incorporated by including 
higher order basis functions for each subsystem boundary 
operator T y . This makes it possible to resolve the fine 
structure of the dynamics and its correlation as well as 
effects due to non-uniform damping over typical scales 
of the subsystem. As one increases the number of basis 
functions, a smooth interpolation from SEA to a full ray- 
tracing treatment is achieved. The maximal number of 
basis functions needed to reach convergence are expected 
to be relatively small thus making the new method more 
efficient than a full ray tracing treatment - in particular in 
the small damping regime. Typical dimensions of T l ° are 
determined by escape-, correlation- and damping-rates of 
the ray dynamics in subsystem j. A priori or a posteriori 
bounds for the size of the basis set needed can thus be 
obtained from dynamical properties of the underlying ray 
flow. 

Representing the ray dynamics in terms of finite di- 
mensional transition matrices may be regarded as a re- 
fined SEA technique. The new method takes advantage 
of the efficiencies of SEA, but includes additional infor- 
mation about the ray dynamics where necessary, thus 
overcoming some of the limitations of SEA and putting 
the underlying SEA assumptions on sound foundations. 
This gives rise to to the name Dynamical Energy Anal- 
ysis (DEA). Note that, like SEA and ray tracing, the 
method is purely based on a classical ray picture and 
is thus inherently a short wavelength approximation, ft 
does not take into account wave like phenomena; from 
a wave asymptotics point of view, these are contained in 
the off-diagonal contributions in Eq. . Wave effects of- 
ten become important in mechanical structures contain- 
ing elements with short and long wave lengths (at the 
same basic frequency) and hybrid SEA - finite element 
methods have been developed in this cas o 19 i 27 ' 30 . An ex- 
tension of these methods to DEA will be of importance. 



V. A NUMERICAL EXAMPLE: A COUPLE TWO-PLATE 
SYSTEMS 

The method has been implemented numerically for a 
coupled two-plate system; the vibrational energy distri- 
bution has been calculated using DEA for plates of dif- 
ferent shape where the coupling between the plates is 
achieved by choosing simply supported boundary condi- 
tions (BC) along a common line of intersection. We as- 
sume clamped BC at the outer edges, that is, Snell's law 
of reflection applies and no losses occur at the bound- 
aries. The two plates have the same thickness and are 
homogeneous otherwise. The BC at the intersection in- 
troduces reflection and transmission and acts as a bar- 
rier thus providing a natural boundary for dividing the 
system into two distinct subsystems. Some of the con- 
figurations considered are shown in the insets of Fig. [5] 
Estimates for the vibrational energy induced by a point 



source in subsystem 1 will be obtain by using DEA and 
will be compared to standard SEA results. 

A. Set-up 

The plates are treated as two-dimensional systems and 
a Fourier basis both in position and momentum space is 
thus an adequate choice for the set of basis functions, 
that is 

^ n ( S ,p s ) = -^ =e 27ri(n 1S /i*+n 2 p,/2) ! 

v 2£i 

with n = (ni,ri2), n%,n2 integers, and s € [0, Li),p s £ 
(—1,1), where is the length of the boundary and i = 1 
or 2. The wave number \p\ cx y/uj is set equal to 1 here. 
Note, that the energy distribution is expected to be fre- 
quency independent as neither the ray paths nor the re- 
flection coefficients at the ray splitting boundary depend 
on k. We also assume for simplicity that the damping 
coefficient a is independent of the driving frequency lu. 
The transmission probability at the intersection of the 
two plates yields for simply supported BC 

Wt (d) = ^cos 2 e 

with 8 e [— 7r/2, 7r/2], the angle between the incoming 
ray and the normal to the surface. 

Given the start and end point s' , s on the boundary 
of either plate 1 or 2, the rays, their lengths and the an- 
gles of intersection (and thus the momentum components 
tangential to the boundary) can be obtained easily, and 
the integral representation of the boundary operator in 
the form (fT7j) will be used. Writing out the Jacobian 
\ds/dp'\, one obtains 

nm J asasw L (s\sj) n m 

where L(s l , s J ) is the length of the trajectory. The weight 
function is given as 

W V = w\ 3 e~ aL 

with a, the damping coefficient, and the reflec- 
tion/transmission coefficients are 

ij , i j\ _ \ $ij ^ s * £ Bj 

w » (s ' s ' ~ \ Sa + (-i)£.u>t(0V, *»')) if s l e B) , 

where B\ denotes the part of the boundary in the coor- 
dinate system s l lying on the intersection of plate 1 and 
2. 



B. Numerical results 

The plates considered in this study all consist of sets of 
straight boundaries^ - such polygonal shapes are typical 
for many engineering applications. Three different set- 
ups have been chosen (see Fig. [5]): 
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FIG. 2. Ratio of mean energy densities in plate 1 and 2 versus damping coefficient a for the three different plate geometries; 
lower right hand corner: relative difference (in percent) between results from SEA (N = 0) and DEA with N = 6 for the three 
configurations. 



• configuration A comprises two sub-systems of ir- 
regular shape with a line of intersection relatively 
small compared to the total length of the bound- 
aries; the two subsystems are thus well separated 
and SEA is expected to work well. 

• configuration B consist of two plates where the 
line of intersection is of the order of the size of 
the system; the only dynamical barrier is posed by 
the BC itself. The standard SEA assumption of 
weak coupling and a quasi-stationary distributions 
in each subsystem may thus be violated. (This con- 
figuration has also been studied i n 26 i 28 ). 

• configuration C has a left-hand plate with reg- 
ular features and rays are channelled out of this 
plate effectively introducing long-range correlations 
in the dynamics thus again violating a typical SEA 
assumption. In addition, the source is chooses at 
the far end of plate 1 in contrast to the other two 
configurations with a source placed close to the in- 
tersection. 

Note, that SEA results are in general insensitive to the 
position of the source, whereas actual trajectory calcu- 
lations may well depend on the exact position especially 
for strong damping and for sources placed close to or far 
away from points of contact between subsections. 



Numerical calculations have been done for finite basis 
sets up to 7i i, ri2 = —N,...N with N < 6. This gives rise 
to matrices of the sizes dimT = 2(2N + l) 2 with basis 
functions covering position and momentum coordinates 
uniformly in both subsystems. Energy distributions have 
been studied as function of the damping rate a. Note, 
that in the limit a — ► 0, the matrix T has an eigenvalue 
one with eigenvector corresponding to an cquidistributed 
energy density over both plates, see the discussion follow- 
ing Eq. (|15p . In the case of no damping, the ray dynamics 
explores the full phase space uniformly on the manifold 
H{X) = uj 2 in the long time limit. Eq. (pt)|) is singular 
for a = and the solutions become independent of the 
source distribution po for a — *■ 0. One obtains 

lim 4~ = lim — = 1 

where pi denotes the mean ray density in plate i averaged 
over the area of the plate and is the corresponding 
mean energy density obtained from Eq. ([3]). 

Results for the relative energy density distribution for 
the two-plate systems are shown in Fig. [21 Increasing the 
basis size - indicated here by the index N - leads to fast 
convergence as is evident from the figures. An SEA-likc 
treatment corresponds to N = 0, here. The lower right 
hand panel in Fig. [5] also shows the difference between 
an SEA and DEA treatment (the latter with N = 6). 
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FIG. 3. Energy density distributions for the three configurations for a = 1 and N = 6; left: wave energy induced by reflections 
from the boundary; right: total wave energy distribution including direct contributions from the source. 



SEA works remarkably well for configuration A, for 
which the main SEA assumptions, that is, irregular 
shape, well separated subsystems and relative small 
damping, are fulfilled. The deviations between SEA and 
the high-resolution result N = 6 are of the order of a few 
percent. Given that SEA describes the energy densities 
here in terms of a system of only two coupled equations, 
this clearly shows the power of SEA compared to, for 
example, ray tracing methods. In configuration B, the 
division into two subsystems is less clear-cut and devia- 
tions from SEA due to the strong coupling between the 
plates may be expected. Indeed, one finds a higher en- 
ergy density in plate 2 than expected from SEA - energy 
dissipates into plate 2 before an equilibrium distribution 
is attained in plate 1. The effect is here of the order of 
10 % and thus still relatively small. 



However, it is not too difficult to devise plate configu- 
rations where significant deviations from SEA occur. In 
conf. C, plate 1 has rectangular shape thus acting as an 
effective channel for transporting wave energy from plate 

1 to plate 2; the plate shape thus induces long range 
correlations and memory effects into the ray dynamics. 
In addition, the source has been placed away from the 
intersection magnifying both the influence of correlation 
effects as well as short range effects due to absorption 
for large a. One indeed finds more wave energy in plate 

2 than expected from an SEA treatment for small a - 
a clear sign that plate 1 acts as an effective channel. 
For large a, however, the wave energy gets damped out 
before reaching the intersection due to the relative long 
path lengths caused by the position of the source. Thus, 
less wave energy than expected from an SEA treatment 
reaches plate 2 and the ratio ei/e 2 is above the SEA 



curve. Note, that the deviations between DEA and SEA 
are now significant and in the region of about 50%. 

Within DEA, it is also possible to resolve the wave 
intensity distribution within each of the subsystems. 
The boundary distribution obtained from Eqs. (|14j) and 
(|20[) can be mapped back into the interior using the op- 
erator Cg 1 in each subsystems, see Eq. (ll"5j) . The spa- 
tial resolution of the wave energy density contains impor- 
tant information about, for example, the (acoustic) radia- 
tion characteristics of sub-elements in the high frequency 
limit. In Fig. [3j typical intensity distributions are shown 
for the three plate configurations at a medium damp- 
ing rate a = 1 and N = 6. The left hand panel shows 
the wave distribution induced by rays reflected from the 
boundary (indirect contributions), the right hand panel 
also includes the direct rays emanating from the source 
point. (It is worth clarifying that only the indirect sig- 
nals have been considered in the results presented in Fig. 
[2]) . The wave intensity plots confirm the observations de- 
scribed earlier; while for configuration A and B, one can 
identify a quasi-equilibrium distribution in each subsys- 
tem characterised by a plateau-structure in each of the 
two plates, the correlated dynamics in configuration C 
leads to a smooth decay of the signal within plate 1. 



VI. CONCLUSIONS 

We have shown that ray tracing methods and SEA are 
closely related and that the latter is indeed an approx- 
imation of the former by smoothing out the details of 
ray dynamics within individual subsystems. We propose 
a numerical technique which interpolates between SEA 
and full ray tracing by resolving the ray dynamics on a 
finer and finer scale. This is achieved by expressing the 
dynamics in terms of linear boundary operators and rep- 
resenting those in terms of a set of basis functions on the 
boundary. The resolution of the dynamics is now deter- 
mined by the number of boundary functions taken into 
account. 

We provide a derivation starting directly from a 
short wave length approximation of the wave equation 
and leading all the way to setting up the basic DEA 
equations; we thus offer a step by step account of 
the approximations and simplifications made. The 
basic SEA assumptions can be tested systematically 
by relating them back to aspects of the ray dynamics. 
Furthermore, extending SEA to DEA allows to enhance 
the range of applicability of an SEA like treatment and 
will lead to robust and numerically efficient tools for 
determining energy density distributions in complex 
mechanical structures. 
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To be precise, polygonal shapes as considered here lead to 
pseudo-integrable dynamics which is strictly speaking not 
even ergodic; at the level of approximation considered here, 
the decay of correlation in the dynamics of irregular poly- 
gons is sufficient in principle to test the SEA assumptions. 



trajectories. If the ray dynamics is chaotic, that is, the 
ray solutions show exponential sensitivity to initial con- 
ditions, one finds that the number of trajectories going 
from ro — > r increases exponentially with their length 1 . 
Regular dynamics such as the solution of the Eq. (|A3|) 
for rectangular or circular geometries leads to a power 
law increase in the number of ray solutions from ro — > r. 



APPENDIX A: RAY DYNAMICS 

A ray or classical dynamics associated with a wave 
equation ^ can be obtained via an Eikonal approxima- 
tion writing the solutions in the form of a phase S(r) and 
amplitude A(r); assuming that the amplitude A changes 
slowly on the scale of the wave length, one obtains a gov- 
erning equations for the phase S alone. For example, for 
the Helmholtz equation with H = c 2 V 2 , one obtains 



2 (V5) 2 



(Al) 



where c denotes the wave velocity (assumed to be con- 
stant here). Dissipative terms are usually incorporated 
in the equation for the amplitude A. The Hamilton- 
Jacobi equation (|A1[) can be solved by the method of 
characteristics. After defining the wave number vector 
p = VS 1 (where we adopt the notation of classical me- 
chanics where p refers to momentum) and the Hamilton 
function 



H(p, r) = c 2 p 2 = lu 2 



(A2) 



one obtains the ray-trajectories (r(r)),p(r)) from Hamil- 
ton's equations 

r = —r = V P H = 2c 2 p; p=-^p = -V r H. (A3) 
dr ' ■'- 



dr 



The fictitious time t is conjugated to the 'energy' uj 2 
and is related to the physical time by t — 2ojt. The 
dimcnsionlcss action S is given as 



S{r,r ) = / dr'p(r') 



(A4) 



where the integration is taken along a ray from ro to r 
on the manifold H(r,p) — u 2 . For homogeneous media 
(c = constant), as considered here, on obtains S = \p\L 
with L(r, ro), the length of the ray path from ro — ► r. 

The ray dynamics in mechanical structures consisting 
of coupled sub-systems will typically entail reflection on 
boundaries, partial reflection/transmission at interfaces 
between two media and multi-component ray dynamics 
including mode conversion. The latter may occur be- 
tween pressure and shear "rays" at boundaries for typi- 
cal boundary conditions (such as free boundaries); note, 
that the different wave components have different local 
wave velocities and will thus follow different equations of 
motion (|A3[) . 

The number of different rays starting in ro (with ar- 
bitrary momentum) and passing through r increases (for 
fixed tu) rapidly with the length or the action of the ray 



APPENDIX B: DERIVATION OF THE RAY-TRACING 
EQ. © 

It will be shown here that Eq. © is equivalent to the 
ray tracing equations ([H} , (fT3"]) . For further details on the 
derivation, see alsoi Starting point is Eq. (flU)) (where 
we set to = 1 here to simplify the notation), that is, 

p(r,r 0> w) = J drj dp' 5(r-<p T r (p' ,r )) 6(u 2 -H(p' ,r )) 

(Bl) 

We write the (^-functions in the form 
5(r-tf(j/, r )) 5(uj 2 -H{ p \ r )) = ]T ±6(^)6^'-^) 

(B2) 



where the index j counts all possible solutions of 

V?(ro,Pj) =r; H(p' j ,r ) = u 2 . 

These are the rays emanating from the source point ro 
and reaching the final point r on the manifold H = to 2 . 
The Jacobian D is 



D = 



d(r,H) 



d{ P ',t) 



dr dH 

dp' dp' 

dr ffH_ 

dt dt 



Making use of the equation of motion (|A3|) . one iden- 
tifies dH/dp' = r' and we have dH/dt = 0. It is 
now convenient to switch to a local coordinate systems 
r = (ru,r±);p — (p\\,p±) at the initial and final point 
where r», pn point along the trajectory in phase space. 
One obtains 



D 



dr 



o 



o 



r r 



dr\ 



dp± 



A (9) 



(B3) 



where A* 9 ' is the geometric contribution to the wave am- 
plitude, Eq. ©. Combining Eqs. {B2]) and (JB3J) with 
(jBlj) . one obtains the diagonal term (jSJ). 



